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Abstract 

This work is a follow-up to the paper, ’’Numerical Relativity as a Tool for Studying the Early 
Universe”. In this article, we determine if gravitational waves can be accurately extracted from a 
dynamical spacetime using an averaging process as opposed to conventional methods of gravita¬ 
tional wave extraction. We calculate the normalized energy density, strain and degree of polariza¬ 
tion of gravitational waves produced by a simulated turbulent plasma similar to what was believed 
to have existed shortly after the electroweak scale. This calculation is completed using two numeri¬ 
cal codes, one which utilizes full General Relativity calculations based on modified BSSN equations 
while the other utilizes a linearized approximation of General Relativity. Our results show that the 
spectrum of gravitational waves calculated from the linear code are nearly indistinguishable from 
those calculated from the nonlinear code using an averaging process. This result validates the use 
of the averaging process for gravitational wave extraction of cosmological systems. 
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I. INTRODUCTION 


Recent work by the Bicep 2 Collaboration |T], has resulted in increased interest in the exis¬ 
tence of primordial gravitational waves. Although these results where later disproven Eli, 
speculation about the characteristics of primordial gravitational waves is still very much 
alive [201 EMU- Our ultimate goal is to determine the characteristics of these waves, given 
different cosmological theories, and if possible, hnd ways to observe them. However, before 
we can do that we must hrst determine how to extract gravitational waves produced by 
turbulence in the relativistic plasma which dominated the universe shortly after the elec- 
troweak phase transition. Can it be accurately calculated using linearilized approximations 
of gravity? Does the chaotic nature of the General Relativity’s nonlinearities signihcantly 
affect the solution isi? Can an averaging process be used to quickly extract gravitational 
waves from a dynamic spacetime na? 

Primordial plasma turbulence is believed to have occurred as a result of stirring caused 
by bubble wall collisions and other chaotic events during inflation and the Electroweak phase 
transition. The characteristic velocity of the turbulent eddies was calculated to be as high 
as 0.65 EH E3]- The magnetic helds around this time may have been as high as 10^° G. 
This number was determined by extrapolating the accepted upper limit on cosmic magnetic 
helds, 10“® G to a much smaller scale factor and assuming that the B oc. relationship 
holds [TH] . 

The author’s previous work [TH| showed how the framework of numerical relativity could 
be used to study the cosmology of the early universe. In this article, we expand on the 
techniques presented and determine the importance of including nonlinear effects. We utilize 
two General Relativistic Magnetohydrodynamic (GRMHD) codes in order to simulate a 
turbulent plasma similar to that of the universe when it was less than a second old. One 
of these codes utilizes full general relativity modeled using a modihed version of the BSSN 
equations in order to simulate the dynamical spacetime background while the other relies 
on a linearized approximation of gravity. We characterize the gravitational waves produced 
and then compare the results produced by both codes. 
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II. OVERVIEW OF THE SOFTWARE 


Our code uses SI units for input and output and geometerized units, c = h = G = are 
used for all calculations within the code. Throughout this article, we will refer to time and 
distance in units of meters. 

As described in the preceding paper [18], the codes used here were specihcally devel¬ 
oped to study early universe cosmology. These codes are based on the Cactus Framework 
(www.cactuscode.org). Cactus was originally developed to perform numerical relativistic 
simulations of colliding black holes but it’s modular design has since allowed it to be used 
for a variety of Physics, Engineering and Computer Science applications. It is currently be¬ 
ing maintained by the Center for Computation and Technology at Louisiana Sate University. 
Cactus codes are composed of a flesh (which provides the framework) and the thorns (which 
provide the physics). The codes used within this work, FixedCosmo and SpecCosmo, are a 
collection of cactus thorns. They are written in a combination of F90, C and C-I--I-. 

Both codes utilize the relativistic MHD evolution equations proposed by Duez [TT] . Both 
codes are also designed to utilize a variety of different differencing schemes including 2nd 
order Finite Differencing, 4th order Finite Differencing and Spectral Methods. The key 
difference between the codes is that while one is capable of solving Einstein’s Equations 
directly (through a modified BSSN formulation) as well as the relativistic MHD equations, 
the other solves the relativistic MHD equations but simulates an expanding spacetime and 
estimates the gravitational wave background without solving directly Einstein’s Equations. 
This allows us to perform a test to determine under what conditions it is important to 
spend computational resources to solve Einstein’s equations directly and if the gravitational 
waves are being correctly extracted from the nonlinear code. The codes were thoroughly 
tested [18] and found to accurately model known GRMHD dynamics. These tests included 
MHD waves induced by gravitational waves test, the consistency of cosmological expansion 
test and shock tests. 

The initial data used was derived from work done by several projects involving primordial 
magnetic fields, phase transitions and early universe cosmology in general [TG] |25l EU ES] 
EH- This study models an extremely high energy epoch of the universe shortly after the 
Electroweak phase transition when the universe was about 10“® s old. The authors chose 
this as the starting point for our study because it was the beginning of the Hadronic Epoch 
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of the early universe. 


III. EVOLUTION EQUATIONS 

The MHD equations used to evolve both numerical codes were based on Duez’s evolution 
equations m. 


dtp* + dj{p*v^) = 0, 

(1) 

dtf + di(a'^y/^ - p*u*) = s. 

(2) 

dtSi + T/) = T°‘^gap,i, 

(3) 

dtB^ + - v^&) = 0. 

(4) 


Here p* is conserved density, is velocity, r is the energy variable. Si is the momentum 
variable, s is the source term, a is the lapse term, 7 is the determinate of the three metric 
and is the stress-energy tensor. The tilde denotes that the term was calculated with 
respect to the conformal metric. 


P* = a^/jpou°, (5) 

f = (6) 

(7) 

s = a+ T°^)dia (8) 


The hrst equation comes from conservation of baryon number, the second derives from 
conservation of energy, the third is conservation of momentum and the fourth is the magnetic 
induction equation. For this simulation we use Geodesic Slicing, a = 1.0, /?* = 0.0 for both 
codes. 

The nonlinear code utilizes a hrst order version of the BSSN equations to simulate the 
background space-time. For hxed gauge conditions, the modihed BSSN equations as dehned 
by Brown [7] are: 
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OqK — O' I Aij + — ] + 47ra(p + S') 


^00 = -^K , 

0 

0o0i = -\oiDiK - K‘^Ci , 
0 

^O'yij 2:0!A^j , 


doA, = 


a{Rij - SnSij) - 2aD{i(t)j) + Aa(pi(pj + AT^ij{2a(pk) 

k 


TF 




do'ykij j^Aij 


kij 


doA^ = --aD^K + 2a (+ GA^^j - Sirf^Sj 
o 


( 9 ) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 

(15) 


The bar denotes a derivative taken with respect to the hducial metric (dehned here to 
have a determinate of one) and the tilde again denotes a derivative taken with respect 
to the conformal metric. Also, C* and Vkij are constraint equations and and are 
proportionality constants, p, S, Sj and Sij are source terms as found in the standard 
version of the BSSN equations. Brown et al also defined: 


Ci = (j)i- DA = 0 , 

Dkij dfcij Dk'^ij 0, 

AT= -A {iklj + llkj — Ijkd , 

+ 7 “[ 2 Af'”K,Aq,„, + Af“,iAf„„] . 

These equations allow gravitational waves to appear organically from the turbulent 
plasma field. For the linear code, we approximate the effect of an expanding spacetime 
and determine the gravitational wave spectrum without utilizing full general relativity. We 
do this by solving the Friedmann Equations and the linearilzed gravitational wave equa¬ 
tions [25] . 


doa = aH 

OqH = -H'^ - + 3p) 
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Where hij are gravitational perturbations calculated by the linearized code and qij are 
their time derivatives. 

IV. CALCULATION OF THE GRAVITATIONAL WAVE SPECTRUM AND PO¬ 
LARIZATION 

The Gravitational Wave Spectrum is output as characteristic strain, energy density and 
degree of polarization by the numerical code [22] • For the full GR code, we utilize an 
averaging process to calculate the gravitational waves produced by the simulation HU. Here 
gravitation perturbations can be found by subtracting the mean value of the metric, gij from 
it’s value at any point and correcting for the scale factor. 


^ij iSij 9ij) /^ 


(16) 


Unfortunately these perturbations are not in the transverse-traceless gauge so the energy 
density must be calculated as 



The brackets denote the average over several wavelengths. By utilizing Geodesic Slicing 
conditions the last two terms are identically zero. The time derivatives of the perturbations 
can be rewritten in terms of the extrinsic curvature, K, lapse, a, Hubble Parameter, H = 
and scale factor, a. 


dohijdoh^^ - ^dohdoh = 2[{2a^KijK^^ + AaHK + QH^) - {aK + 3Hf]/a^ (18) 


The mass density is then normalized by dividing by the critical density 


3H^ 


(19) 


8irG' 


This results in a normalized energy density. 


6 




Qn = -{{2a^KijK^^ + AaHK + 6H^ - {aK + (20) 

The normalized energy density for the linearized code is then simply calculated as 


rir — 




( 21 ) 


For both codes, the magnitude of the characteristic strain is calculated using the quadratic 
mean of the plus and cross polarizations of the gravitational waves perceived to be traveling 
along the x, y and z axes. The plus polarizations are calculated by finding the difference 
between the diagonal transverse perturbation terms while the cross terms relate directly to 
the perturbations in the off-diagonal transverse terms. 


^zz) 

(22) 

hy '^{hzz hxx) 

(23) 

hz — ~{hxx ~ ^yy) 

(24) 


(25) 

— h 

'^y ’^XZ 

(26) 

hz = hxy 

(27) 


(28) 

/ hx{k, ty + hy{k, ty + hz{k, ty 

(29) 


According to Kahniashvili 133, the degree of polarization can be dehned as, 

{h^*{k,t)h^{k',t) + h^*{k,t)h^{k' ,t)) 

Where the right and left polarizations for waves traveling along the x, y or z axis are 
defined as. 


hf = hf + ih. 


K = K- 
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(31) 

(32) 






By expanding the left and right polarizations into plus and cross polarizations, the po¬ 
larization degree can be rewritten as, 




2{Im[hl{ki, t)]Re[h^ {ki,t)] - -Re[h+ {ki, t)]Im[hf {ki, t)]) 

{\ht{h,t)\^ + \hUh,tW) 


(33) 


V. EXPERIMENTAL SET-UP 

Our simulation begins shortly after the electroweak energy scale, at the beginning of the 
Hadronic Epoch. The numerical values for the initial conditions where calculated using 
the available literature [2S] and are the exact same for both codes. These calculations 
give us an initial temperature of 1.30 x 10^^ K at time 1.0 x 10“® s. The scale factor 
and Hubble Parameter are a = 2.096 x 10“^^ and H = 7.99 x 10®s“^ respectively. The 
critical mass/energy density at the time was 1.14 x 10^^^. The characteristic velocity of 
the turbulent eddies for all runs was set to 0.25 [3T1 [33]. The magnetic held at the time 
could have been as large as 10^^ G so we set the amplitude of the magnetic helds to 10^^ G 
for our simulations. 

We ran three sets of simulations using the nonlinear and linear codes with the exact same 
initial conditions. The hrst simulation utilized a grid 1000 m cubed with 64^ grid points 
and a timestep of 10“® m. These dimensions were chosen so that any resulting gravita¬ 
tional waves would correspond to the frequency range of Pulsar Timing observations once 
universal expansion is taken into account. The other two used higher or lower resolutions in 
order to establish convergence in our results. Geodesic slicing conditions, periodic boundary 
conditions and Fourier spectral differencing were used for all simulation runs. There were 
no shocks or discontinuities in the system so we did not utilize our HRSG routines. We also 
used a 3rd order Iterative Grank Nicolson time scheme for time integration. All runs started 
with no initial gravitational waves but the density, temperature, velocity and magnetic helds 
were all perturbed to model that of an early universe space-time |aiini[2Il[221ESlEH|. The 
initial density and temperature were perturbed by random phase cosine functions with an 
amplitude proportional to their wavenumber, k, ehectively a Fourier series. The initial 
magnetic held consisted of random phase cosine waves with an amplitude proportional to 
their wavenumber squared, The initial velocity held consisted of random phase cosine 
waves with an amplitude proportional to their wavenumber cubed, k^. Each run utilized 64 



TABLE I. Degree of Polarization at the final time. 


X - Direction Y - Direction Z - Direction 

Nonlinear (Mean) 

0.0 

-1.56125E-17 

0.0 

Linear (Mean) 

0.0 

0.0 

0.0 

Nonlinear (Standard Deviation) 

0.592164099 

0.518579627 

0.569916566 

Linear (Standard Deviation) 

0.592333947 

0.515687302 

0.572444906 


processing cores on the Maxwell computing cluster at the University of Houston’s Center 
for Advanced Computing and Data Systems. Over 10,000 iterations were produced for each 
data run. 


VI. RESULTS 

It should be noted that these results have not been extrapolated for present day obser¬ 
vations and only represent a relatively short period in the early universe. The strain and 
normalized energy density outputs for all runs appeared to be independent of k, so we chose 
to focus on the mean value of these quantities. As one can see from Figures 1 and 3, the 
difference between the average strain as calculated by the nonlinear and linear codes is neg¬ 
ligible. We see that the same is true for the Energy Density as shown in Figure 2. Based on 
the data, we can assume that strain and energy density calculations derived from the linear 
code is accurate to within a part in a thousand of those derived from the nonlinear code. 

As shown in Table 1, The degree of polarization for the gravitational waves appears to 
be the same for both the linear and nonlinear calculations. This calculation was performed 
from a sample of 5 different times between 1.0005 x 10“®s and 1.0006 x 10“®s. As expected, 
since the initial data was symmetric, there does not appear to be any bias in the degrees of 
polarizations for the X, Y and Z directions or between left and right polarized gravitational 
waves generated by the turbulent system. Further, the statistics of the polarizations for 
both the nonlinear and linear simulations appear to be the same to two signihcant hgures. 
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Strain vs Universe Age 


1.2E-12 



Normalized Gravitational Wave Energy Density vs Universe Age 


l.OOOOE-06 l.OOOlE-06 1.0002E-06 1.0003E-06 1.0004E-06 1.0005E-06 1.0006E-06 



—^GRMHD 
-SRMHD 


Universe Age 


FIG. 1. Strain of gravitational waves produced by the linear and nonlinear codes. Normalized 
Energy Density of gravitational waves produced by the linear and nonlinear codes. 


VII. DISCUSSION 

In this article, we tested the averaging process to extract gravitational waves from nonlin¬ 
ear cosmological simulations. We did this by using cosmological linear and nonlinear codes 
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%diff vs Universe Age 
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FIG. 2. The %diff was calculated by dividing the difference in strain (energy density) by the strain 
(energy density) of the nonlinear code. 


to solve the exact same problem and comparing the results. We looked at three character¬ 
istics of the resulting gravitational waves, strain, normalized energy density and degree of 
polarization. We found that the results agreed for all three measures. The tiny differences 
that did occur may be partially explained by small nonlinear effects. We believe that these 
results prove that the averaging process outlined in Section 4 of this paper is an accurate 
method of extracting gravitational waves for cosmological systems that lack spacetime sin¬ 
gularities. This holds true for gravitational waves with strains up to 10“^^ and normalized 
energy densities as high as 10“®. Further work is needed to test the limits of the averaging 
process. Also, as a result of this work, we can conclude that the linear and nonlinear codes 
produce the same gravitational waves to within a part in 10^. 
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